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Abstract 

Computing the viability kernel is key in providing guarantees of safety and proving existence of safety- 
preserving controllers for constrained dynamical systems. Current numerical techniques that approximate 
this construct suffer from a complexity that is exponential in the dimension of the state. We study 
conditions under which a linear time-invariant (LTI) system can be suitably decomposed into lower- 
dimensional subsystems so as to admit a conservative computation of the viability kernel in a decentralized 
fashion in subspaces. We then present an isomorphism that imposes these desired conditions, particularly 
on two-time-scale systems. Decentralized computations are performed in the transformed coordinates, 
yielding a conservative approximation of the viability kernel in the original state space. Significant 
reduction of complexity can be achieved, allowing the previously inapplicable tools to be employed for 
treatment of higher-dimensional systems. We show the results on two examples including a 6D system. 



1 Introduction 

Constrained dynamical systems have received a tremendous amount of attention due to the presence of safety 
constraints and hard bounds that appear in many practical scenarios. Providing guarantees of constraint 
satisfaction and facilitating synthesis of constraint-satisfying controllers therefore is highly desirable, partic- 
ularly in safety-critical applications. A class of safety-critical systems known as envelope protection problems 
is concerned with ensuring that the trajectories remain in a safe, bounded "envelope" (subset) of the state 
space for a given time horizon. Such problems arise in e.g. flight management systems [l]-[4] where the safety 
constraints are defined as the aircraft's aerodynamic envelope and consequently the system must ensure 
that certain combinations of states are avoided to prevent stalling or other undesirable behaviors. Other 
application domains include control of depth of anesthesia 5 , aircraft autolanders 6 , automated highway 
systems control of under-actuated underwater vehicles [8^, stockout prevention of storage systems in 
manufacturing processe s [o] , and management of a marine renewable resource 10 , to name a few. 



Viability theory 11-13 provides a set- valued perspective on the behavior of the trajectories inside a given 



set. Thus it is naturally suited to handle envelope protection problems. By duality, minimal reachability 



14 is also capable of analyzing such problems by investigating the behavior of the trajectories outside of the 
envelope. For simplicity, in this paper we only focus on the constructs generated within the framework of 
viability theory. The viability kernel is the set of initial states for which there exists at least one trajectory of 
the input-constrained system that respects the state constraint for all time. It is shown in [12] and (by duality 
in y^) that the viability kernel is the only construct that can be used to prove safety /viability of the system 
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and to synthesize inputs that preserve this safety; cf. [Tg", Chap. 1-2] for more detail. In general an exact 
computation of the viability kernel is extremely difficult if not impossible. Instead, approximations of this 
set are computed. Such computations have historically been subject to Bellman's "curse of dimensionality" 
The numerical algorithms that approximate the viability kernel and its associated control laws (e.g., 



17 



14 



18 -20'), collectively referred to as Eulerian methods 15 , rely on gridding the state space and therefore 
their computational complexity increases exponentially with the dimension of the state. This renders them 
impractical for systems of dimension higher than three or four. 

This paper presents a part of our efforts to address the curse of dimensionality by enabling the use 
of Eulerian algorithms for higher-dimensional LTI systems (and by extension, hybrid systems with LTI 
dynamics). We decompose the structure of the system, applying Eulerian algorithms on each individual 
lower-dimensional subsystem in a decentralized fashion. Significant computational gains can be obtained, 
since instead of one costly centralized computation on the full-order system, multiple less expensive sub- 
system computations are performed. The results are then mapped back to the full-order space to obtain a 
conservative approximation (i.e. an under- approximation) of the viability kernel. The contribution of this 
paper is twofold: 1) We investigate various structures on system matrices that must be satisfied so that 
the behavior of the constrained system for envelope protection problems (with simply-connected, compact 
constraints) can be inferred conservatively from subspace decentralized analyses (Section |3|. 2) We then 
present an isomorphism through which the desired structure is imposed on the system (albeit under cer- 
tain conditions) to facilitate decentralized computations in the transformed space (Section [4|. Numerical 
examples are provided in Section [5j 

1.1 Related Work 

Complexity reduction for viability and minimal reachability has been addressed by many researchers. A pro- 
jection scheme in 21 based on Hamilton- Jacobi (HJ) partial differential equations (PDEs) over- approximates 
the projection of the true minimal reachable tube in lower dimensional subspaces, with the unmodeled dimen- 



sions treated as a disturbance. Similarly, 22 decomposes a full-order nonlinear system into either disjoint or 
overlapping subsystems and solves multiple HJ PDEs in lower dimensions. More recently, a mixed implicit- 
explicit HJ is presented in 23 for nonlinear systems whose state vector contains states that are integrators 
of other states. The complexity of this new formulation is linear in the number of integrator states, while 
still exponential in the dimension of the rest of the states. These techniques assume that the system itself 
presents a certain structure that can be exploited. 



In 24 , an approximate dynamic programming technique is presented that, although still grid-based, 
enables a more efficient computation of the viability kernel. The viability kernel (similarly to 25 ) is expressed 
as the zero sublevel set of the value function of the corresponding optimal control problem. It is assumed that 
the value function, which is a viscosity solution of a HJB PDE, is differentiable everywhere on the constraint 
set. The PDE is then discretized and the resulting value function is numerically computed on a grid using 
a function approximator such as the /c-nearest neighbor algorithm. The error-bounded approximation is 
not conservative (it is an over-approximation) but converges to the true viability kernel in the limit as the 
number of grid points goes to infinity. 

Another related approach is the search for a barrier certificate 26 , a Lyapunov-like function that forms 
a separating hyper-surface between any two given sets A and B in the state space. If there exists a function 
non-positive on A and positive on S, and whose Lie derivative (along the vector field) is non-positive on its 
zero level set for all states and controls, then no trajectories will ever go from A to B. This technique can 
be adapted to analytically describe the boundary of the infinite-horizon viability kernel: A certificate must 
now be formulated such that at every state along its zero level set there exists a control that makes the Lie 
derivative non-positive. For systems with polynomial vector fields and semi-algebraic constraints, efficient 
techniques based on Sum of Squares can be used to find the barrier certificate^ 

■"^This method cannot be used to formulate the finite-horizon viabihty kernel which may be useful when, for example, the 
infinite-horizon kernel is empty, or when safety is to be verified/enforced over a finite time interval. Moreover, there are no 
guarantees that a barrier certificate can be found for a given system no matter how simple its dynamics (even when a Lyapunov 
function is already known). 
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Recently, we presented a connection between the viability kernel and efficiently-computable classes of 
reachability constructs known as maximal reachable sets. Owing to this connection, scalable numerical 



algorithms (collectively referred to as Lagrangian methods [15 ) such as 27 33], originally developed for 
maximal reachability, can now be used to approximate the viability kernel. We presented two algorithms for 
LTI systems with convex constraints based on piecewise ellipsoidal representations [5^ and support vectors 



34 that have polynomial complexity. In contrast to these results, the technique presented here reduces 
the complexity indirectly by decentralizing computations. The benefit of this approach is that it allows 
useful features of Eulerian methods such gradient-based control synthesis and handling of arbitrarily shaped 
nonconvex constraints be taken advantage of. 



2 Problem Statement 

Consider the continuous-time system 



x = f{x,u) (1) 



with state space A' := (a finite-dimensional vector space), state vector x{t) G A', and input u{t) G U 
where U is di compact (closed and bounded) and convex subset of R^. The vector field f : x U ^ is 
assumed to be Lipschitz in x and continuous in u. Let 

^[o,t] •= {u: [0,t] W measurable, u{s) e U a.e. 5 G [0,t]} . (2) 

With an arbitrary, finite time horizon r > 0, for every t G [0,r], xq G A', and u{-) G ^[Q^t]^ there exists a 
unique trajectory (.xq.u'- [0^^] ^ ^ that satisfies ([T]) and the initial condition (,xq,u{^) = ^o- 

For a nonempty, simply-connected, compact state constraint set /C C A' we are concerned with computing 
the following backward construct}^ 

Definition 1 (Viability Kernel). The finite-horizon viability kerne^ of K is the set of initial states for which 
there exists an input such that the trajectories emanating from those states remain in JC for all time t G [0, r]: 

Viab[o,,](/C,ZY) := {xo G /C | 3u{-) G Vt G [0,r], ^,„,(t) G JC} . 

Initial states belonging to this set are viable under ([T]), and the corresponding control laws are safety- 
preserving. The powerful Eulerian methods are capable of directly computing the viability kernel and 
its safety-preserving control policies. However, they rely on gridding the state space, and therefore are 
computationally intensive. Although versatile in terms of ability to handle various types of dynamics and 
constraints, the applicability of these techniques has been historically limited to systems of low dimensionality 
(up to 4D in practice) due to their exponential complexity. 

We restrict ourselves to LTI systems of the form 

x = Ax^Bu (3) 

described by the matrix notation 

S:=[A\B] (4) 
with constant, appropriately sized A and B matrices. 

Problem 1 (Decentralized Viability), i) Identify a structure on A and B for which the viability kernel can 
be conservatively reconstructed from its subsystem analyses, ii) Find an isomorphic state space for ([3| in 
which the system has this desired structure. 



^By duality the arguments presented in this paper also hold for the minimal reachable tube of /C^; cf. 
"^The infinite-horizon viability kernel Viab]^+ (/C, is also known as the maximal controlled- invariant set l35|. 
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2.1 Preliminaries 

Notation For a set ^ C A', and 2-^ denote the complement and the power set of ^ in A', respectively. 
For brevity, ||-|| denotes the infinity norm. For a constant matrix A = [aij] G W^^'^ the induced norm is 
||A|| := sup^^^r. ,y_^o ~ ^^'^i<j<n'^^i\(^ij\' For a Lebesgue measurable function / : R ^ defined 
over an interval [ta^h] we denote ||/|| := ^f,] = ^^^te[ta,tb] 11/(^)11 < ^ linear transformation 

of 5 in Q using a nonsingular matrix T G R^><^ is defined as 5' = T-i(5) := [ T'^AT \ T'^B ] . A linear 
transformation of a set ^ C A' under the same mapping is 3^ = T~^A {y \ y = T~^a^ a G A}. 



Definition 2 (Disjoint Input). The input u = [ui ■ ■ ■ Up]^ (ZW is disjoint across two subsystems 



xi = AiXi + A12X2 + Biu, 
±2 = A2X2 + A21X1 + B2U 



of an LTI system with xi G R^ and X2 G 



du. 



^ ifyse{l,...,p},iy^j, 
0, 



dBiU , ^ dBjU 
7^ ^ 



dus 



(5a) 
(5b) 



(6) 



andU =Ui XU2, where Ui is any (possibly degenerate) subset ofW 
acting directly on subsystem i draws its values. 

Definition 3 (Unidirectionally Coupled). The subsystems 

xi = AiXi + Biu, 

±2 = A2X2 + A21X1 + B2U 



from which the portion of the vector u 



(7a) 
(7b) 



with disjoint input across them are said to be unidirectionally coupled since the trajectories of (7b) are 



affected by those of (7a)^ while (7a) evolves independently from (7b). The worst-case unidirectional coupling 



can be characterized by ||A2i||. 

Definition 4 (ETUC). A subsystem is said to be externally trivially uncontrollable (ETUC) if it possesses 
a null input matrix. 

Remark 1. The condition onU in Definition\^ enures that the inputs acting on each subsystems are inde- 
pendent of one another. This condition is satisfied for most physical systems where actuators are commonly 
uncorrected, or for a system with an ETUC subsystem (in which case the shape ofU becomes irrelevant). 
In the most general case, however, U can be (under-) approximated by a cross-product set. 



3 Decentralized Viability Computation 

We begin by arriving at the desired structure on system matrices that would allow for decentralized (and 
conservative) computation of the viability kernel. Throughout the paper we assume a partitioning of (|4| 
that results in two subsystems. The arguments can be easily generalized to multiple subsystems as discussed 
in Section [421 



3.1 Why Decoupling of A Alone is Insufficient 

Consider the following system with block diagonal A-matrix, and a 5-matrix of generic form: 

ueU. 



Xi 




'Ai 




Xi 






_X2_ 




A2_ 




_X2_ 







Denote the two subspaces of K" in which the subsystems evolve as 

Si := R'^ and §2 := K""''- 



(8) 



(9) 
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Let UiX be the projection of the vector x = [xi X2]^ ^ ^ onto S^: 

Uix = Xie Si, (10) 

and UiJC the projection of the set /C C A' onto S^: 

UiJC = {xi eE>i\3xeJC, UiX = . (11) 

Lemma 1. For any t and u{-) G ^[o,^] the projection of trajectory ^ of system ([8| with initial condition 
^xo,n(0) = Xq is a subsystcm trajectory initiating from the projection of xq: 

n^UAt) = eu,.oJt)- (12) 

Proof. Hi [il] = Hi ([^0^ I] [^l]u) = AiUiXi + BiU. □ 

Corollary 1. 

U.u(t) elC^ ^k..o,u{t) e n,/C. (13) 
Later we wih show and utiUze the fact that under certain conditions this impUcation is bidirectional. 

Proposition 1 (Wrong Approximation). For dynamics (|8| the cross-product of subsystem viability kernels 
of projections of JC is a superset of the viability kernel of K: 

Viab[o,r](/C,ZY) C Viab[o,r](ni/C,ZY) x Viab[o,r] (n2/C,ZY). (14) 

Proof. 

xo e Viab[o,r](/C,ZY) ^ 3u{-), Vt, ^xoA^) ^ ^ 

Uixo e Viab[o,r](ni/C,^/) A U2X0 G Viab[o,r](n2/C,^/) 
^ xo e Viab[o,r](ni/C,^/) x Viab[o,r](n2/C,^/). 

□ 

The fohowing counter example demonstrates that an inclusion in the opposite direction does not hold for 
system ([8|; That is, Viab[o,r](^5^) 2 Viab[o,r](ni/C,^/) x Viab[o,r](n2/C,^). Consider the point = [}] and 
constraint set JC = [—1, 1] x [—1, 1]. We seek to compute the viability kernel of this set under the dynamics 
Xi = xi -\- u and X2 = ^2 — u and input constraint u G [—1, 1]. The point x' belongs to the cross-product 
of subsystem viability kernels (since subsystem 1 can use u = —1 while subsystem 2 can use u = at the 
same point to keep Hix' in H^/C), but does not belong to the actual full-order kernel (since no input exists 
that can keep the system in JC). As such, when the system is in the form of ([8| performing the analysis 
on subsystems would yield an over- approximation of the viability kernel. This stems from the fact that the 
input is non-disjoint across the subsystems. On the other hand, we do have the following correct inclusion 
even with a non-disjoint input. 

Lemma 2. The following holds for system 

Viab[o,,](/C,ZY) D (Viab[o,,]((ni/C^)^ZY) x S2) U (Si x Viab[o,,]((^2/C^)^ZY)) . (15) 
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Proof. 



xo e (Viab[o,,]((ni/C^)^ZY) X §2) U (Si X Viab[o,,]((^2/C^)^ZY)) 

=^ (V^i(-), G Ei/C^ A eu^.o^t) ^ n2/C^))' 

^ Xo e Viab[o,r](^5^)- 

Definition 5 (Ill-Posedness). We say that a viability problem is ill-posed if the state constraint is empty. 



□ 



Proposition 2 (111- Posed Approximation). When K is a bounded subset of X (which is the case in most 
envelope protection problems) the approximation in Lemma^is ill-posed. 

The proof should be clear from the fact that for any bounded set JC we have {HiJC^y = 0. 
3.2 Suitable Structures for Decomposition 

Consider a system with block- diagonal A-matrix and a 5-matrix that ensures a disjoint input across the 
subsystems, for instance 

xi] \Ai 0] \xi] \Bi 0] \ui] . . 

X2J [ A2J [X2\ [ B2\ [U2\ ^ ' 

when U =Ui XU2. 

Assumption 1. The set JC is a cross-product of two (arbitrarily- shaped) sets in S^. 

Corollary 2. Under Assumption^ the projection of a trajectory is contained in a set if and only if the 
subsystem trajectories are contained in the projection of the set: 

U«We/c^4,„,„^(t)eni/c. (17) 

Theorem 1. The viability kernel of JC under ([I6]) can be computed exactly using subsystem kernels: 

Viab[o,^](/C,ZY) = Viab[o,^](ni/C,ZYi) x Viab[o,^] (02/0,^^2). (18) 

Proof. 

xq G Viab[o,r](^7^) ^ Vt, ^ {t) eJC 

^ 3u{-), Vt, {^h.xoA^) ^ ^ ^kxoA^) ^ n2/C) (by Assumption (TJ 

^ 3u,{-), Vt, eu^.o,uAt) ^ ni/C A 3u2{-)^ Vt, ^l,,,^u,{t) ^ n2/C (via disjoint input) 
^ Uixo G Viab[o,r](ni/C,^/i) A 02^0 G Viab[o,r] (n2/C, ^/2) 
<^ Xo G Viab[o,r](ni/C,^/i) x Viab[o,r] (112/0,^/2). 



□ 
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Remark 2. The use of any decomposition technique for correct (conservative) approximation of the viability 
kernel is contingent on satisfaction of Assumption [7] as shown previously. When JC does not satisfy this 
assumption, it can be under- approximated by the union of direct- product sets. The viability kernel can be 
computed for each set separately in lower dimensions (which increases the computational complexity only 
linearly in the number of sets). The union of the resulting kernels in full dimensions under- approximates 
the true viability kernel. Parallelization of viability calculations in each subspace could further reduce the 
computational time. 

In general, we may not be able to simultaneously obtain a decoupled A-matrix and a disjoint input. 
Instead, suppose that the system is of the form 



Xi 




'A, 







Xi 






_X2_ 




A 


A2_ 




_X2_ 








(19) 



which automatically ensures that the input u is disjoint across the subsystems regardless of the shape of 
U since one of the two (unidirectionally coupled) subsystems is ETUC (Remark [T]). This system can be 
rewritten as 



Xi 




'Ai 




Xi 











_X2_ 




A2_ 




X2 









A 



(20) 



The evolution of xi is completely independent of the evolution of X2. Its effect on the lower subsystem, 
mapped through A, can be viewed as an exogenous input to the lower subsystem, that takes values on the 
(possibly time- varying) subset V(-) of the upper subspace Si. Treating this additional input in the worst-case 
fashion results in conservatism. Hence, define the following construct: 

Definition 6 (Discriminating Kernel). Consider a system with adversarial inputs: control u{t) G U and 
disturbance v{t) G V{t), where V: [0,r] 2^^"" is a point-wise convex and compact set-valued map from [0,r] 
toRP-. Let 

^[o,t] '= [0,^] measurable, v{s) e V{s) a.e. s e [0,t]}. 

To be conservative, we assume non-anticipative strategies p for one of the inputs^ The finite-horizon dis- 
criminating kernel of JC is the set of initial states for which there exists a control such that the trajectories 
emanating from those states remain in JC for every disturbance for all time t G [0, r] : 

Disc[o,,](/C,ZY, V(-)) := {xo G JC \ 3p: 1^o,r] ^ ^[o,r], V^(-) G 1^o,r], Vt G [0,r], UAv]At) ^ 



We will use a subscript to distinguish a construct formed under (20) when xi for the lower subsystem 
is treated as an adversarial disturbance. 



Lemma 3. The viability kernel of a set JC under ( |2Q[ ) is a superset of the discriminating kernel of JC when 
xi is treated as a worst-case disturbance (assumed to draw values from some time-varying set V(-) point-wise 
convex and compact in Si) to the lower subsystem: 

Viab[o,,](/C,ZY) D Disc[o,,](/C,ZY, V(-))*. (21) 

Proof. Let ^ denote the trajectory of the system when xi is treated as a disturbance to the lower subsystem. 



Xo G Disc[o,r](^,^, V(-))* ^ 3pM(-), Vv(-), Vt, 4o,pM,^W ^ ^ 

^3ix(-),Vt, e,,,,(t)G/C 
=^ Xo e Viab[o,r](^7^)- 



(a specific disturbance) 



□ 



ma p p : '^o,t] ~^ '^[o,t] is non-anticipative for u if for every ^'(•), v' {-) G ^[o,t]7 ^{^) — ^\^) implies = a.e. 

s G [0, t] Note that for linear systems the Isaac's condition holds [m], and therefore it does not matter which input is 

selected to play with non-anticipative policies. 
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Definition 7 (Invariance Kernel). Consider a system with a disturbance input v{t) G V{t) as its only input, 
where V(-) is defined as in Definition^ The finite-horizon invariance kernel of a set JC is the set of initial 
states that remain in IC for every disturbance for all time t G [0, r]: 

Inv[o,,](/C, V(-)) := {xo & IC \ M') e y[o,r], Vt e [0,r], UA*) & -^l- 



Theorem 2 (Main Decentralization Result). The viability kernel of a set JC under ([T9| can be conservatively 
approximated using the subsystem viability /invariance kernels as 

Viab[o,,](/C,ZY) D Viab[o,,](ni/C,ZY) x Inv[o,,] (Hs/C, V(-))*, 

where V: [0,r] ^2^^; t^ Viab[o,r-t](ni/C,ZY). (22) 

Proof We first show that the inclusion holds for any set P C Si in which xi takes value. Since both inputs 
(control u and "disturbance" v := Xi e V) are disjoint across the two subsystems we have 

Viab[o,,](ni/C,ZY) X Inv[o,,](n2/C,P)* = Disc[o,,](ni/C,ZY, {0})* X Disc[o,,] (Hs/C, {0}, P)* 

^^Disc[o,,](/C,ZY,P)*. (23) 

With V = V(-), inclusion ( [22| ) follows from Lemma [sj 

Viab[o,,](ni/C,ZY) X Inv[o,,](n2/C, V(-))* = Disc[o,,](/C,ZY, V(-))* C Viab[o,,](/C,ZY). 

Note that the set- valued map V(-) at time t is the finite- horizon viability kernel of the upper subsystem 
over the interval [0,r — t]. This map is continuous (it is both lower and upper semicontinuous (cf. [12]) at 
every point in its domain) and non-decreasing 19 (i.e. V{t) 3 V(s) Vt G [5,r], 5 G [0,r]), with UiJC being 
its upper-limit in the sense of Kuratowski (Definition^ as t ^ r~ and Viab[o,r] (ni/C,^/) its lower-limit as 
t 0+. Furthermore, since IIi/C and U are convex ana compact and the dynamics linear, the sets V{t) are 
also convex and compact at every t. From this we have that Inv [o, r-s] (n2/C, V(5)) is continuous, convex and 
compact for every 5, and non-decreasing over s G [0,r] 12 . 

We use these statements to argue that a digression from the formulation in (22) loses its sufficiency to 
guarantee an under- approximation in the sense that if the uncertainty set is assumed to be a subset of V{t) 
for any t then the cross-product may not generate an under- approximation of the viability kernel: Consider 
a set-valued map V(-) s.t. 3i G [0,r], V(£) C V{i) (e.g. a constant set Viab[o,r] (ni/C,^/) Vt). It is clear from 
([23| that 

Viab[o,,](ni/C,ZY) X Inv[o,,](n2/C, V(-))* 5 Disc[o,,](/C,ZY, V(-))* (24) 

since for any set C, Inv^Q (C, V(£))* C Inv^Q (C, V(f))* and therefore Inv[o, r-s] (C, V(-))* C Inv[o, r-s] (C, V(-))* 
Vs G [0,f]. There is no guarantee that this superset in ([24| is a subset of Viab[o,r](^7^); Lemma [s] is no 
longer applicable. On the ffip side, if V(-) is such that V{t) 3 V{t) for any t G [0,r] (e.g. a constant set HiJC 
Vt), then an excessively conservative under-approximation is obtained. □ 



3.3 Sub-Interval Formulation and Decentralized Algorithm 

In practice, we can perform the analysis over sub-intervals (similarly to [37]) while still maintaining con- 
servatism. During each sub-interval the set V(-) is sampled and kept constant in backward time. Such 
sub-interval analysis is possible via the semi-group property in both subspaces as well as the following 
results in S2. 

Proposition 3. For N := r/q, N time steps each of length q G we have that 

N-l 

f] CInv[o,,](n2/C,V(-))* (25) 

i=0 

where Ci = Inv[o,g] (C^+i, V((i + 1)^))* with Cn = II2/C. 
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Proof. Notice that since {V{t)}^^Q is a non-decreasing sequence of compact and convex sets with V{t) C 
Si =: R^^ we have that for a fixed for every i, V((i + l)q) 5 V{t) \/t G [0, {i + l)g]. Using this, the fact 
that Ci C Ci+i ^Cn and the semi-group property we have 

nAT-l 
d^Mie [0, N -I], \/vi{-) e {vi : [0, q] measurable, 
i=0 

v,{s) e V{{i + l)q) a.e. s G [0,g]}, Vt G [0,g], G C,+i 

^ Vz G [0, - 1], Vvi(-) G {^i : [iq, {i + 1)(7] M^" measurable, 

G V{s) a.e. 5 G [ig, {i + l)g]}, Vt G [iq, {i + G C^+i 

^ Vv(-) G {v: [0,r] ^ R^^ measurable, v{t) G V(t) a.e.}, Vt G [0,r], ^^^^^(t) G Cat 
^ xo G Inv[o,r](n2/C, V(-))*, 

where is the concatenation of functions Vi over [0,r]. 



□ 



In the limit this set converges to the invariance kernel with unsampled input set. 



Definition 8 (Kuratowski upper and lower limits 19 ). Let {A{s)}ses be a sequence of subsets in a metric 
space (E^d). The upper-limit of A{s) as s ^ s is 



Lim sup A{s) := < X ^ E \ lim inf d{x^ A{s)) = ^ , 

s^s L s^s 

where d{x^A) := inf^^^ a). Its lower-limit is 

Lim inf ^(5) := \ x e E \ lim d{x , A{s)) = ' 



Proposition 4. Denote by Cn{q) CliLo^ intersection of N = r/q sub-interval invariance kernels 

from Proposition^ For the sequence of subsets {Cn{q)}q>o we have 

LimsupCn(<7) = Inv[o,r] (n2/C, V(-))*. (26) 

Proof. Given q^ define a piecewise constant set-valued map Vsh(^;<?) •= V(z(7) Vt for which i is the unique 
integer in {1, ... , N} satisfying t G {{i — l)g, iq] when t varies backwards from r to (i.e. a backward sample 
and hold of V(-)). Recall that V(-) is non-decreasing and continuous, and V{t) compact for every t. Clearly, 
Vsh(-; q) 5 V(-) Wq. The sequence {Vsh(-; Q)}q>o converges to V(-) from outside: We say that v{-; q) G Vsh(-; q) 
iff v(t',q) G Vsh(t;^) Vt. As ^ ^ 0+, Vi;(-;g) G Vsh(-;^) Ve > Vt B{v{t',q),e) H V(t) 7^ 0, where B(x, e) 
denotes the ball (associated with a metric d) of radius e centered at x. In other words, \/v{-]q) G Vshi'^Q)^ 
3v{-) G V(-) s.t. limsupg^o+ d{v{-),v{-; q)) = liminf<^^o+ d{v{-),v{-; q)) = 0. So limg^o+ d{v{-),Vsh{-] q)) = 0, 
and therefore Liminfg^o+ Vsh(-;^) = ^(•)- On the other hand, we know from the semi-group property that 
Cn{q) = Inv[o,r](n2/C,Vsh(-; <?))*• Hence, 

LimsupCn(<?) = LimsupInv[o,r](n2/C, Vsh(-; <?))* = Inv[o,r] (n2/C, Liminf Vsh(-; <?))* = Inv[o,r] (n2/C, V(-))*. 

□ 

Using this formulation we can perform the decentralized analysis in Theorem [2] via Algorithm [T] over 
sub-intervals. 
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Algorithm 1 Sub-Interval Decentralized Computations 



1: N ^ r / q > Assumed integer. 

2: Cn ^ ns/C 

3: Vat ^ Hi/C 

4: for i = TV - 1 to do 

5: Ci ^ Inv[o,g](Ci+i, Vi+i)* 

6: Vi ^ Viab[o,g](Vi+i,ZY) 

7: end for 

8: return Vo x Co > ^ Viab[o,r](^7^) 



3.4 Bounding the Approximation in §2 

Notice from Theorem [2] that the computed construct in the upper subspace is exact in that 

Hi Viab[o,r] (^5 ^) = Viab[o,r] (ni/C, U) . 



(27) 



On the other hand additional conservatism is introduced in the lower subspace S2 due to treating the effect of 
the upper subsystem as a worst-case disturbance. Quantifying this error remains an open problem. However, 
we can formulate a qualitative lower bound on the shrinkage of the invariance kernel in S2 in backward time. 
This bound will be expressed in terms of system-specific (and ultimately, design-specific) parameters that 
form the desired structure (19): 



Following 1 38 , the invariance kernel in S2 can be expressed as 



Inv[o,,](n2/C,V(-))* = n (e-'^'^2JCe f e-'^^'AVit 



r)dr 



te[o,r 



(28) 



with denoting the Pontryagin difference. Let B{S) be the norm-ball of radius S G about the origin, 
and define 7^: R+ ^ R+, 

e^ll^2|| _ 1 

vis) := ^j^. (29) 

Bounding the contribution of the uncertainty (disturbance) in computation of the invariance kernel over the 
interval [0,^] we have 37 that 



Jo 



-rA- 



AVie - r)dr C B 



(ir 



-rA- 



AV{e - r)dr 



e'-||A2ll ii^i 



sup ||a;||(ir 

xeV{9-r) I 



Cb(\\A\\ sup 11x11 fe^ 

\ xevie) Jo 

CS(||A|| sup \\x\\rj{0)] 
Clearly, this contribution is weakened as ||A|| ^0. Further, we have 



(30) 
(31) 
(32) 
(33) 



N-1 



N-1 



n n e-^^^C,+ieS ||A|| sup \\x\\rj{q)] C Q Inv[o,,] (C,+i, V((z + l)g))* (34) 



=0 \te[0,q] 



xeviii+i)q) 



with Cn '= n2/C. From the dual of the results in 137 , we know that the Hausdorff distance of the two sets 



in the inclusion above decreases as g ^ 0+, and tends to zero if V{iq) = S(sup^ 



eviiq) 



11x11). The Kuratowski 
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upper-limit of the left-hand-side of (34) is therefore Inv[o,r] (n2A^, as g ^ 0^ (via Proposition [4|. Now, 



notice that for sufficiently small <7 <C 1, 



2 



r]{q) = lim ^ ^" '"^ < Jim ^ ^ ^ ^ 7, ^ = ^ + ^2) + 0{q'), 35 

where ^(742) and n = dim(S2) respectively denote the largest singular value and the dimension of the 
lower subsystem. Therefore provides a qualitative lower-bound on how much Inv[o,r] (n2/C, V(-))* can 
shrink in backward time in terms of n, the magnitude of the unidirectional coupling ||A||, the supremum 
of V{t) (the viability kernel in Si), and the largest singular value ^(^2) of the lower subsystem. If we can 
choose h appropriately, assign the slow eigenvalues to the lower subsystem, and weaken the effect of the 
disturbance (uncertainty) as much as possible by minimizing ||A||, we can expect the conservatism to be 
reduced considerably. The proposed modified Riccati transformation in Section [4] provides this fiexibility 



while imposing the desired structure (19) on the system. 



3.5 Decentralized Viability in Transformed Coordinates 

Suppose that for a general system ([3| under which a centralized viability computation is known to be 
burdensome, there exists an invertible transformation z = T~^x such that in the new coordinates the system 
S = T~^{S) has the form of either ([l6| or (19). Suppose that Assumption [l] is satisfied for T~^JC. When 



the transformation yields decoupled ^d-matrix as well as disjoint input. Theorem [T] under the transformed 
dynamics S becomes: 

Corollary 3. Viab[o,r](/C,ZY) = T Viab|^^](T-i/C,ZY) = T (viab|^^] (niT-i/C,ZYi) x Viab|^^](n2T-i/C,ZY2)); 
where the superscript S is used to specify when a construct is formed under the transformed dynamics. 

For the more general case Theorem [2] implies: 

Corollary 4. Viab[o,^] (/C, ZY) 3 T (viab|^](niT-i/C,ZY) x Inv|^^] (n2T-i/C, V(-))*) withV{t) := Viab|^^_,](niT-i/C,ZY) 
Vt e [0,r]. 

Decentralized analysis over sub-intervals are performed similarly to Algorithm [l] and a lower-bound for 
the shrinkage of the invariance kernel in S2 can be formulated according to ( [34| ) with Cn = 112 T~^/C and the 

respective transformed system matrices. Note that in Si, EiT"^ Viab[o,r](^7^) = Viab^^](niT~^/C,Z//), 
and that the computed construct in S2 is a guaranteed under- approximation of the projection of the actual 

viability kernel in that subspace, i.e. Inv^^] (n2T~^/C, V(-))* C n2T~^ Viab[o,r](^7^)- We present one such 
transformation next. 



4 The Riccati-Based Transformation 

We draw upon the so-called Riccati transformation — a two-stage coordinate transformation based on the 
solutions of a nonsymmetric algebraic Riccati equation (NARE) and a Sylvester equation. This transforma- 
tion, originally introduced in ^39j for decoupling of singularly perturbed systems, was later generalized in |40] 
to larger classes of autonomous LTI systems. An in-depth overview of the application of this transformation 
in optimal control theory, singular perturbation theory, and asymptotic approximation theory can be found 
in ^41j, while more recent advances are given in [42^j43j . 
Let (|4| be partitioned as 

^ ^ r ^11 A12 
^21 ^22 



B2 



(36) 
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with All e ]R^><^, Ai2 e r^x(^-^), A21 e r(^-^)x^, A22 e ]R(^-^)x(n-/c)^ ^ ^kxp^ ^ 
for some k < n. Now consider the nonsingular transformation matrices 

T2 

where In denotes the n x n identity matrix. With L G m("^-^)x^ and M G M^x(n-fc) ^YiqX satisfy 



h 





-L 


In — k 


h 


M ' 





^n—k 



(NARE:) 
(Sylvester:) 

the transformed system is 

S' = T-^{S) 



^{L) := LAn - A22L - LA12L + A21 = 0, 
y{M) := (An - Ai2L)M - M{A22 + LA12) 



il2 



0, 



All -A12L 

All - A12L 




^12 
A22 + LA12 

A22 + LA12 



Bi 
LBi + ^2 



(/ - ML)Bi - MB2 
LBi + B2 



-k) xp 

(37) 
(38) 



(39) 
(40) 



(41) 
(42) 



Solutions to (39) and (40) may not always exist. The above procedure is referred to as the (standard) Riccati 
transformation. If the control input is disjoint across the subsystems of S" (and thus the transformation 
imposes a structure similar to (16)), Corollary [s] can be employed to approximate the viability kernel in a 
decentralized fashion based on subsystem analysis. 

4.1 The Modified Riccati Transformation 

For the more general case, on the other hand, we propose the following transformation that imposes a 
structure given in (19) which also relaxes the condition on the shape of the set U. Corollary [4] can thus be 
employed to compute a conservative approximation of the true viability kernel. 

4.1.1 Transformation 1 (ETUC Subsystem) 



Consider a transformation through which the lower subsystem can be made ETUC. That is, in (41) for the 
transformation matrix Ti we seek an L in ^(L) that is also a solution of LBi ^ B2 = 0. 



Assumption 2. "^(^J) C ^{B^), where ^{X) is the column-space of matrix X . 

). Under Assumption^ the class of solutions of LBi = —B2 w.r.t. L G m(^-^)x^ can be 

C:= \^-B2B\^ Z - ZBiBI, ZgM^^"^)^'^! 



44 



45 



Lemma 4 ( 

characterized by 

with t denoting the Moore-Penrose pseudoinverse. 



(43) 



Assumption [2] is the necessary and sufficient condition for solvability of LBi = — i?2- Substituting (43) 
for L in ^(i) we obtain 



^{Z) := + r + z(Ai2 - BiB\Ai^Z{BiB\ - I) 

+ {A22 - B2BIAi2)z{BiBI - I), 



(44) 
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where 



E = -{BiBl - I)(^An+Ai2B2Biy 

r = (A22B2-B| + - B2BI (^Ai2B2BI + An). 



(45) 
(46) 



To eliminate the non-invertible term {BiB\ — I) from the right-hand side of (44) we equate M(Z) to 
some rank correcting term 8J^{Z) with 

^{Z) := z[Ai2 - BiB\Ai2)z+ [A22 - B2B\Ai2)z (47) 

and 5 e K\{— 1,0} a finite (but possibly large) parameter such that {BiB\ — {6 + 1)1) is nonsingular: 

^{Z) = Z5 + r + z(Ai2 - SiSj A12) Z{BiB\ - I) 

+ (A22 - B2B\Ai2)z{BiB\ - I) 



= Z5 + r + ^{Z){BiB\ -i) = s^{z). 



(48) 
(49) 



Simple algebraic manipulation and post-multiplication of ^{Z) — 5^{Z) = by (^i^j — ((^ + 1)/) ^ results 
in a NARE in the variable Z: 



^i(Z) := Zin - A22Z - ZA12Z + ^21 = 



(50) 



with ill = S {BiB\ -{5^ 1)1) \ A21 = r {BiB\ - ((5 + 1)/) \ A12 = {BiBIAi2 - A12), and A22 = 
{B2BIA12-A22). 

Proposition 5. If a root Z G ^i^-^)^^ of the NARE (50) exists, it constitutes an L e C that simultaneously 
satisfies 



LBi + ^2 = 0, 

^{L) = LAii - A22L - LA12L + A21 = S^{Z). 



(51a) 
(51b) 

□ 



Proof. By virtue of (49), a matrix Z that satisfies ( |50[ ) also satisfies ([51| via (43). 
Remark 3. If p > k the set C reduces to the singleton {— ^2^1} ctnd the method still applies. 

Theorem 3. The transformation (37) with L G m(^-^)x^ obtained through Proposition^ makes the lower 
subsystem in (36) ETUC. Moreover, the coupling terms are altered such that the effect of the upper subsystem 
on the evolution of the lower subsystem is parameterized by 5. 

Proof. 



S' = T-\S) = 



All - A12L A12 
LAii - A22L - LA12L + A21 A22 + LA12 



Bi 
LBi + B2 



All - A12L A12 
S^{Z) A22 + LA12 



Bi 




(52) 
(53) 
□ 



Remark 4. Note that the imposed 5 -parameterization of the off-diagonal term S^{Z) in (53) provides 
an additional degree of freedom in adjusting (minimizing) the coupling of the two subsystems in the new 
coordinates. This will be discussed further in Section \4.1.3\ 
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Nonsymmetric Riccati equations have long been an active area of research . To solve (50) we draw 
on the fixed-point algorithm described in 40 and derive the necessary conditions for the existence and 
uniqueness of a real root Z. Suppose (52^1^12 — ^22) is invertible. Define initial values as 



:= (BaSl A12 - A22) 'r{B^Bl-{5 + l)l) \ 
Ao ■.= E{B,Bl-{5+l)iy' -{B,bIAu-Au)Zo. 



To find Z we look for 
by solving 



D:= Z-Zo 



^i{D) := DAo - (b2B\Ai2-A22 + Zq{B^bIAi2 - Ai2))d 

- D{BiBIAi2 - Ai2)D + ZoAo = 0. 
Lemma 5 ( 40, Lem. 1]). Suppose (^2^1^12 — ^22) '^s nonsingular. If 

\\{B2B\Ai2-A22y^\\< ^ 



(54) 
(55) 

(56) 
(57) 



then (57) has a unique real root D that satisfies 

0<\\D\\ < 



3(||Ao|| + ||BiBI^i2-^i2||||^o||) 



2Po||||^o| 



(58) 



BiB\Ai2-Ai2 



\Z4 



(59) 



and is the fixed-point solution of the contraction Dk+i = ^i{Dk) given by 
^i(Z)fc) := {B2BIA12 - A22y' (^ZoAo + DkAo 



Zo{BiBIAu - Au)Dk - Dk{BiB\Ai2 - ^12)^^)- 



(60) 



Remark 5. As in it can he shown that the relative error '= \\Df. — D\\ / \\D\\ after k iterations is 
bounded above by 



ek< (^\\{B2B\A^2- A22) '||(Po|| + ||SiBMi2-Ai2||||Zo||) 



(61) 



and decreases as \5\ increases since \\Aq\\ and \\Zq\\ are inversely related to \5\. 

For a given ^, using Dq = as initial condition we compute D iteratively. The fixed-point solution 
= ^1(1^*) is then used to obtain Z = + Zq which in turn solves ^i{Z) = in (50) and results in a 
matrix L, through (43), that satisfies both equations in (51). 

4.1.2 Transformation 2 (Unidirectionally Coupled Subsystems) 

Consider the NARE 

^2{M) = {An - Ai2L)M - M{A22 + LA12) - M{S^{Z))M + A12 = 0. (62) 
For a given L, (5, and Z, if there exists a solution M that satisfies ([62|), we obtain the following: 



Theorem 4. The transformation (38) with M eR^^^^ satisfying NARE (62) makes the subsystems ir 
(53) unidirectionally coupled. 
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Proof. 



All - A12L - M8^{Z) M^^k^^ 



Bi 




(63) 



□ 



Remark 6. In the transformed coordinates the lower subsystem remains ETUC. Furthermore, the 5 -parameterization 
of the unidirectional coupling between subsystems is also preserved. 

Before further analyzing the unidirectional coupling term let us derive the necessary conditions 

for the existence and uniqueness of a solution M to ([62]) to be used with the same convergent iterative 
procedure described previously. For a given ^, Z, and L, let (An — A12L) be invertible and the initial values 
be defined as 



We seek M by forming 
and solving 



Mo := - [All -A12L) ^Ai2, 
No := A22 + LA12 + 6^{Z)Mo. 

J := M -Mo 



h{J) := JNo - [All - A12L - 5Mo^{Z)) J + 5J^[Z)J + MqA^o = 0. 



Lemma 6 ( 40, Lem. 1]). Suppose [An — Ai2L^ is nonsingular. If 

1 



(A11-A12L) 



< 



then (67) has a unique real root J that satisfies 

o< lUII < 



\{\\No\\ + \\S^{Z)\\ \\Mo\\) 



2||iVo||||Mo|| 



(64) 
(65) 

(66) 
(67) 

(68) 



||iVo|| + ||,5^(Z)||||Mo|| 
and is the fixed-point solution of the contraction Jk+i = ^2{Jk) given by 

^2{ Jk) ■■= {An - AuL)'^ (MoTVo + JkNo + 5Mo^{Z)Jk + SJk^{Z)Jk 
Remark 7. The relative error e^ := \\Jk — J\\ / after k iterations is bounded above by 

k 



Bk < 3||(A 



111 



A12L 



^)"'||(l|A^o|| + ||<5^(^)||||Mo||) 



(69) 



(70) 



(71) 



and decreases as \\8^[Z)\\, IIA22II; and ^[Aii—Ai2L) "^|| decrease. This occurs when the ill- conditioning of 
the A-matrix increases (e.g. in the case of two-time- scale systems; see (^7/ and the references therein) and 5 
is chosen such that \\5^[Z)\\ is minimized. 

Using Jo = as initial condition we compute J iteratively. The fixed-point solution J* = ^2(^*) is then 
used to obtain M = J* + Mq which in turn solves ^2(M) = in (|62). 

Note that both conditions (ISSl) and ( 68 ) are conservative and their satisfaction ensures rapid convergence 



(usually within 2 or 3 iterations). In practice, the right-hand-side of these inequalities can be relaxed up to 
10 times in most cases without causing divergence. 
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4.1.3 The Unidirectional Coupling Term (Choosing S) 

Finally, we analyze the unidirectional coupling term S^{Z) and its behavior with respect to the free pa- 
rameter 5. Since Z is an implicit function of we adopt the extended notation S^{Z{S)) to reflect this 
dependency. First, we formalize a conservative upper-bound on || as an explicit function of S. 

This assures that the unidirectional coupling remains bounded for almost all admissible values of the free 
parameter 6. 

Proposition 6. The worst-case unidirectional coupling between the two subsystems in the transformed co- 



ordinates^ i.e. \\5^{Z{5))\\ in (63); is (conservatively) bounded above such that 



^ Mlkr^J ^nkr^J'' v.eR\{-i,o}, (72) 

where the constants a and b are independent of S and are determined by a := a{b/P)'^, b := 3||5i5| ||7/3; 
7:= ||r|| 11(^22 -B2B|^i2)"'||, a := \\Ai2-BiBlAi2\\, and (3 := P22 - 52-81^12 1|. 

Proof. The proof is provided in the Appendix. □ 



Now consider inequalities (58) and (68), which are dependant on S. Adequately chosen and sufficiently 
large values of 5 help ensure that these conditions are met. On the other hand, choosing S exceedingly large 
defeats the purpose of ^-parameterization of the unidirectional coupling term, since it can be shown that as 
6 grows, ||^^(Z((5))|| approaches a problem-dependant constant that may not necessarily be an extremum 
point. 



Proposition 7. lim \\S^{Z{S))\\ = ||r|| with T given by (46). 

(5^±oo 



Proof. This proof is also provided in the Appendix. □ 

It follows from Proposition [7| that < Ms \\S^{Z{S))\\ < \\T\\. Therefore naively letting \S\ oo 
essentially removes the added flexibility associated with the ^-parameterization in the modified Riccati 
approach and instead enforces a trivial solution L = —B2B\. While for some systems this solution may yield 
the smallest possible unidirectional coupling between the resulting subsystems (i.e. a unidirectional coupling 
with the least infinity norm), in most cases a carefully chosen 5 not only facilitates the satisfaction of the 



convergence conditions rt58k and (68), but also further minimizes the worst-case unidirectional coupling. 



Thus, formulated as an optimization problem, we seek a 8 that solves the following: 

subj. to ([58| and (|68|. 

Note that this is a nonconvex problem, and in general, / may be a non-smooth function of S. However, a 
global optimum need not be computed. Any suboptimal solution can be used as long as that solution yields 
a satisfactory degree of unidirectional coupling between the subsystems in the transformed coordinates. An 
approximation to the optimum point can be obtained numerically, for example by fine-griding the real line 
or using the bisection algorithm. 

In practice, while the exact shape of the function / is problem-dependant, we have found (but not proven) 
that in most cases it exhibits a behavior similar to that of an absolute value proper rational function (over 
a discontinuous domain) of the form 

/(5) = |^+ci|+C2, ySeV, (73) 
where V C M\ { — 1,0} is the union of the two segments of the real line for which the magnitude of 5 is large 



enough such that ([58|) and (68) are both satisfied, /c G N, /c : odd, cq = — ci((^*) , 5* = arg min^^;^ /(J), 



02 = minsey f{S), and ci = (lim^^ioo f{^)) - C2 = - C2. 
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Figure 1: The worst-case unidirectional coupling f{S) = ||(5^(Z((5))|| (x's) and its approximation f{S) = 
l^^p^ + 0.55| + 1.82 (dashed) computed for Example [l] The interval (-15, +15) over which ([58| and ([68| 
are violated is labeled as "infeasible region". The asymptote lim^^^ioo /(^) = ||r|| (dash-dotted) is also 
shown. The minimum of f{6) occurs when S ~ +50. 



Example 1. Consider the system 



1.5072 


3.3984 


0.1300 


-0.0884" 




"-0.7433" 


5.0644 


-2.6683 


0.0227 


0.1689 




-2.2528 


0.1156 


-0.1863 


0.5686 


0.2648 


, B = 


-0.9075 


-0.0808 


0.0229 


0.4915 


0.5949 




0.6036 



Fi^. [7] ^/loi^;^ f{S) and its approximation f{S) 
A randomized, empirical test in 



27.65 



-0.55| + 1.82 evaluated where (58) and (68) hold. 
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Section 4.4.2] examines the potential affect of the system dimension 
n on the magnitude of the unidirectional coupling and the amount of time consumed by the decomposition 
process. While the test shows an increasing trend in average values, there is significant variance. In addition, 
the time required for the decomposition process (even for the highest dimension n = 16 in our test) is still 
negligible (~1.5s) compared to the time required for the actual viability computations. 



4.2 Recursive Decomposition 

A recursive decomposition when the standard Riccati transformation can be used is straightforward. Suppose 
that the modified Riccati transformation is used throughout the process. In deeper level recursions, the 
decomposition can be applied to the uppermost subsystem since that subsystem is controlled whereas every 
other subsystem is ETUC. For example, to decompose a 6D system into three 2D subsystems, in the first 
recursion level, the partitioning can be chosen such that the resulting upper (controlled) subsystem is 4D and 
the lower (ETUC) subsystem is 2D. In the second recursion level, if the solutions exist, the 4D subsystem 
is then decomposed into two 2D subsystems. Note that in the recursive application of the decomposition, 
when the modified Riccati transformation is employed, all subsystems but one are ETUC. Therefore, this 
iterated decomposition may result in an excessively conservative under-approximation of the true viability 
kernel. 



4.3 Riccati-Based Viability in Lower Dimensions 

In the new coordinates z = T~^x, T = T1T2, the subsystem dynamics are governed by 

ii = {All - A12L - SM^{Z))zi + Biu, (74) 
Z2 = {A22 + LA12 + S^{Z)M)z2 + B2U + 6^{Z)zi (75) 

with S^{Z) = when the standard Riccati transformation yields disjoint input, or ^2 = when the modified 
Riccati transformation is employed. In the latter case, (5 = (5* is precomputed so as to minimize \\S^{Z)\\. 
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In addition, the transformation automatically assigns the slowest eigenvalues to the lower subsystem. These 
in turn prevent excessive conservatism in approximation of the construct in S2. Analysis over sub-intervals 
are performed according to Algorithm [ij and a qualitative lower-bound for the shrinkage of the invariance 
kernel in S2 can be formulated according to (34) with Cn = Il2T~^JC and A = S^{Z). 



5 Numerical Examples 



Among Eulerian methods we use the Level Set Toolbox (LS) v. 1.1 48 for our analysis. All computations are 
performed on a dual core Intel-based machine with 2.8 GHz CPU, 6 MB L2 cache and 3 GB RAM running 
single-threaded 32-bit Matlab 7.5. 



5.1 4D Cart with Two Inverted Pendulums 



Consider the linearized model of a cart with two separately mounted inverted pendulums from 49, Ex. 2.2.1] 
with h = 30, h = 35: 






1 





0" 







0.3920 





-0.0327 





, B = 


-0.0033 











1 





0.0560 





0.2753 







-0.0005 



The state vector x G consists of angular displacement of each inverted pendulum from vertical and the 
corresponding angular velocities; The input G M, < 10, arises from a force applied to the cart. 

Note that despite the sparsity of the system no permutation matrix can recover our desired structures 



(16) or ([T9| (the graph representation of this system is a strongly connected digraph). We decompose this 
system using the presented Riccati-based technique into two 2D subsystems, with unidirectional coupling 
determined by the solution L = —B2BI regardless of the value of 5: 





0.3920 





0.9524 



0.1429 

0.2800 






1.0500 











-0.0033 











We choose JC such that in the transformed coordinates we have the constraint set JCz := {z \ \\z\\ < 
0.5, z = T~^x, X G JC}. We seek to identify the set of initial states for which there exists a bounded 
control law that keeps the angular displacement of the pendulums contained in JCz and thus within a ball 
of finite radius about their upright positions, despite control saturation. We perform the analysis over 50 
sub-intervals. LS v. 1.1 only accepts hyper-rectangular input sets. To comply with this limitation we modify 
Step [5] in Algorithm [1] so that Ci ^ Inv[o,g] (C^+i, Box(Vi+i))*, where Box(^) is the interval hull of A. 
Conservatism in Proposition |3] is preserved since Box(V(zg')) ^ V(iq). Computations are performed over a 
grid with 41 nodes in each dimension using a first-order accuracy for r = 3s (Fig.[2|. The computation time 
for the actual and the transformation-based kernels were 1098.48 s and 4.27 s, respectively. The Riccati-based 
kernel covers 74% of the volume of the full-order set (calculated based on the number of grids contained in 
each set). 



5.2 Arbitrary 6D System 

Consider the two-time-scale system x = [f^^^ \x^{f^^\u with e = 0.1, and A G R^><^ and 5 G 
matrices randomly drawn from a normal distribution A/'(0, 1): 



A = 



--0.3557 


-0.3078 


-0.6097 


2.0275 


-1.3636 


-0.4131- 




- 1.0720 


-0.8153- 


0.1233 


-1.6441 


0.2404 


-0.6431 


0.0517 


-0.1454 




-1.7390 


-0.7181 


1.8857 


-1.1748 


-1.2502 


-0.7252 


-0.7801 


-0.3972 




-0.8292 


-0.4906 


-0.0194 


-0.0779 


-0.0208 


0.0160 


-0.0465 


0.0535 


, B = 


0.0156 


0.0540 


-0.0486 


-0.0192 


0.0781 


0.1017 


0.0838 


-0.0518 




-0.0960 


0.0875 


_ 0.0043 


-0.0849 


-0.0228 


-0.0901 


-0.0319 


-0.1143. 




.-0.0347 


-0.0054_ 
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Figure 2: Riccati-based (solid, dark) vs. actual (transparent, light) viability kernels in the transformed 
coordinate space for Example 5.1 



We decompose this system into two 3D subsystems using the modified Riccati transformation with 
^* ^ -25: 

■ 1.0720 -0.8153- 
-1.7390 -0.7181 
-0.8292 -0.4906 


. _ 

The constraint /C is chosen such that this set in the new coordinates is a nonconvex set formed by the 
cross-product of the union of a sphere and a hyper-rectangle as shown in Fig. [Sj We choose U such that 
—0.5 < ui < 0.5 and 0.5 < U2 < 1. (The shape of U need not be rectangular since one of the subsystems 
is ETUC.) Decentralized approximation of Viab[o,2](^7^) are carried out over 50 sub-intervals using 151 
nodes in each dimension and a second-order accuracy (Fig. [3|. The overall computation time was Ih 
(including calculation of J*, transformation matrices, the decomposition, and projections which took only 
a few seconds). In contrast, the actual kernel is prohibitively computationally expensive to compute with 
LS for any meaningful grid resolution. Moreover, on average 350 MB of RAM was used in the Riccati-based 
viability calculations (of which 110MB was to store the grid), whereas the computation of the full-order 
kernel would require about 380 TB (terabyte) merely to store the grid. 

5.3 Comparison With Schur-Based Decomposition ( |50| ) 

In [50] we presented a Schur-based decomposition technique that is applicable to almost any LTI system. 
In contrast, the decomposition method presented here is based on two nonsymmetric algebraic Riccati 
equations. The existence of solutions to these algebraic equations, however, is limited by a number of 
conditions on system matrices and is therefore heavily problem dependent. Indeed, as pointed out earlier, 
the conditions are more likely to be satisfied as the ill-conditioning of the original system matrices increases — 
e.g., for two-time-scale systems However, when the algebraic Riccati equations do converge, the resulting 
subsystems could potentially yield less conservative kernel approximations than in the case of the Schur-based 

^cf. [l6| Figure 4.6] for the fraction of tests on randomly generated systems for which a solution existed. 



A" 



--0.3472 


-0.1553 


-0.5243 


0.1252 


-1.6394 


0.2499 


1.8832 


-0.9445 


-1.1162 


0.0069 


-0.1476 


-0.0544 


-0.0523 


-0.0749 


-0.0097 


-0.0015 


-0.0604 


-0.0238 







-0.1011 
0.1474 
-0.1425 







0.0244 
0.0156 
0.0200 



■ 





0.1152 
-0.0571 
-0.0762. 
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Figure 3: The constraint set (transparent) and its Riccati-based viability kernel in 3D subspaces of the 
transformed coordinates for Example [5] 



decomposition; See 16, Section 4.5.4] for a simple example. In general, however, it is the problem under 
study that determines which decomposition method is more suitable. A better strategy may be to use both 
decomposition techniques if possible and take the union of their resulting sets to obtain a more accurate 
under- approximation of the viability kernel than what could be achieved using each individual technique. 



6 Conclusions and Future Work 

We considered the problem of guaranteed safety and constraint satisfaction in moderately-dimensioned, 
safety-critical LTI systems with compact, simply-connected state constraints. To provide such guarantees 
the computation of the viability kernel is required. Historically, the algorithms that approximate this set — 
known as Eulerian methods — are based on gridding the state space. While powerful and versatile, their 
computational complexity increases exponentially with the dimension of the state which renders them im- 
practical for systems of dimensions higher than three or four. We investigated conditions under which the 
viability kernel can be conservatively approximated in a decentralized fashion in lower-dimensional sub- 
spaces. We then presented a new similarity transformation that imposes such conditions on the system, 
thereby allowing us to employ Eulerian methods on higher-dimensional systems. The transformation is best 
suited to two-time-scale systems. 

It is possible (although uncommon) that the transformation matrix can become poorly-conditioned due 
to pseudoinverses and numerical algorithms involved, resulting in the state constraint set in the transformed 
coordinates becoming too severely distorted under the linear map to be of any practical use. An upper- 
bound on the condition number in terms of the system matrices and the free parameter S is provided in 
p^, Appendix B.2]. We are currently investigating possible remedies that would ensure a well-conditioned 
transformation matrix. 



With the particular system structure (19) considered in this paper, the computations in the upper 
subspace are exact. On the other hand, the lower subspace computations are subject to accuracy loss since 
the formulated disturbance is assumed to play optimally at all times, aiming to shrink the construct in that 
subspace. While this is to ensure that we obtain a conservative approximation, in reality it is quite likely 
that the input is not always adversarial. Moreover, here we have only required the disturbance signal be 
measurable, and thus it can vary discontinuously. We know, however, that the trajectories of the upper 
subsystem are continuous in time. Restricting the disturbance input to draw from the subclass of continuous 
signals may result in a more accurate approximation in the lower subspace. In either case, quantifying the 
accuracy loss in Lemma [3] is an open problem. Another future direction is in investigating alternative system 
structures to the ones considered in Section [321 
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Appendix 



Proof of Proposition\S[ From the matrix inversion lemma, (Y+UCV)-^ = Y-^-Y-'^U{C-'^+VY-^U)-'^VY-^, 
with Y = -{S + 1)1, [/ = Bi, C = /, and V = we have 



{B^Bl - (6 + 1)1)-' = -^(/+ is.Bt). 

Using this, (47), ([54|, ( [56| ), ([59|, multiplicative and triangular inequalities, and > 1, 

¥^(Zim < 1*1 WIIZoll + + f!(||Z„|| + IIBID) 



2IW 



oil II^OI 



||Ao||+a||Zo||; ■ -V"""" ' ||Ao||+a||^o| 
<|,5|(9a||Zof + 3/3||Zo||) 

< \d\ (9a7l (B^bI -{d + 1)/) ir + 3/^7|| (^i^l - (<5 + 1) J) 



< 1(51(^907^ 

1 /H^ 



1 



< 



\s\\\s + i\ 



5 + 1 

a - 



i<^i+i 



•3/^7 



1 



+ 



b, V(5 G 



(5+1 
1,0}. 



(76) 



□ 



Proof of Proposition^^ Notice from (58) and (61 ) that for large values of 5, Z can be closely approximated 
by its initial value Zq . Using ( 76 ) , 



hm \\5^{Z{6))\\ 



lim 



^P2Qi(/ + = ||o + P2gi|| = ||r|| 



with Qi := (52^1^12 - A22) 'r. Pi := (A12 - PiPIAi2), P2 := (52^1^12 - A22). 



□ 
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